*
*		PROGRAM dynsimladderplot
*	
*		3/20/16
*		Andy Philips
*		Texas A&M University
* -------------------------------------------------------------------------
* -------------------------------------------------------------------------
* NOTES: --program is identical to paneldynsimpieinter.ado but creates
*		   a ladderplot.
*		 --interaction is only suitable for dichotomous variables.
*		 --interaction variable is interacted with ALL variables (including
*		   dummy variables)
*		 --if dummyset not given, dummy will be set to 0
* -------------------------------------------------------------------------


capture program drop dynsimladderplot
capture program define dynsimladderplot , rclass
syntax [varlist] [if] [in], [ dvs(varlist max = 12) shockvar(varname) 	  ///
Time(numlist integer > 1) SHock(numlist)] 								  ///
[shockvar2(varname) shock2(numlist)] [shockvar3(varname) shock3(numlist)] ///
[shockvar4(varname) shock4(numlist)] [shockvar5(varname) shock5(numlist)] ///
[dummy(varlist)] [dummyset(numlist)] [sig(numlist integer < 100)]		  ///
[range(numlist integer > 1)] [saving(string)] [INTERaction(varname)]	  ///
[INTYpe(string)] [id(varname)] [basetime(numlist)] [endtime(numlist)]

version 8
marksample touse
preserve

if "`sig'" != ""	{						// getting the CI's signif
	local signif `sig'
}
else	{
	local signif 95
}
local sigl = (100-`signif')/2
local sigu = 100-((100-`signif')/2)

if "`range'" != "" {						// How far to simulate?
	local range `range'
}
else {
	local range 20
}
if `time' >= `range' {
	di in r _n "The range of simulation must be longer than the shock time"
	exit 198
}
if "`dvs'" == ""	{
	di in r _n "Must specify dependent compositional variables in dvs( )"
	exit 198
}
if "`shockvar'" == ""	{
	di in r _n "A shockvar, not included in [varlist], must be specified"
	exit 198
}
if "`intype'" == "" & "`interaction'" != "" {
	di in r _n "Intype must either be (on) or (off) in intype"
	exit 198
}
* ------------------------Generating Variables & Run Model ---------------------
local lvars 
local dvars
local LIvars
local DIvars
qui foreach var of varlist `varlist'  {		// get d. and l. indep vars
	tempvar L`var' D`var'
	bysort `id': gen `L`var'' = l.`var'
	bysort `id': gen `D`var'' = d.`var' 
	local lvars `"`lvars' `L`var''"'
	local dvars `"`dvars' `D`var''"'
	if "`interaction'" != ""	{				// and the interacted ones, if needed
		tempvar LI`var' DI`var'
		bysort `id': gen `LI`var'' = l.`var'*`interaction'
		bysort `id': gen `DI`var'' = d.`var'*`interaction'
		local LIvars `"`LIvars' `LI`var''"'
		local DIvars `"`DIvars' `DI`var''"'
	}
}

tempvar lagshockvar diffshockvar			// same for the shock variable
qui bysort `id': gen `lagshockvar' = l.`shockvar'
qui bysort `id': gen `diffshockvar' = d.`shockvar'
local lagshockvariables `"`lagshockvar'"'
local diffshockvariables `"`diffshockvar'"'
qui if "`interaction'" != ""	{					// if the interaction is on
	tempvar LIshockvar DIshockvar
	bysort `id': gen `LIshockvar' = l.`shockvar'*`interaction'
	bysort `id': gen `DIshockvar' = d.`shockvar'*`interaction'
	local LIshockvariables `"`LIshockvar'"'
	local DIshockvariables `"`DIshockvar'"'
}
qui forv i = 2/5	{							// and the additional shocks
	if "`shockvar`i''" != "" {
		tempvar lagshockvar`i' diffshockvar`i'
		bysort `id': gen `lagshockvar`i'' = l.`shockvar`i''
		bysort `id': gen `diffshockvar`i'' = d.`shockvar`i''
		local lagshockvariables `"`lagshockvariables' `lagshockvar`i''"'
		local diffshockvariables `"`diffshockvariables' `diffshockvar`i''"'
		if "`interaction'" != ""	{
			tempvar LIshockvar`i' DIshockvar`i'
			bysort `id': gen `LIshockvar`i'' = l.`shockvar`i''*`interaction'
			bysort `id': gen `DIshockvar`i'' = d.`shockvar`i''*`interaction'
			local LIshockvariables `"`LIshockvariables' `LIshockvar`i''"'
			local DIshockvariables `"`DIshockvariables' `DIshockvar`i''"'
		}
	}
}

local dumIs											// for interacted dummies
if "`interaction'" != "" & "`dummy'" != ""	{
	foreach var in `dummy'	{
		tempvar dumI`var'
		bysort `id': gen `dumI`var'' = `var'*`interaction'
		loc dumIs `"`dumIs' `dumI`var''"'
	}
}
else	{
	loc dumIs ""
}

local model
local i 1
qui foreach var of varlist `dvs' {
	tempvar ldepvar`i' ddepvar`i' mdepvar`i'
	bysort `id': gen `ldepvar`i'' = l.`var'				// gen lag DV
	bysort `id': gen `ddepvar`i'' = d.`var'				// gen diff DV
	bysort `id': gen `mdepvar`i'' = `var'				// grab means
	if "`interaction'" != ""	{			// create interacted lag DV
		tempvar LIdepvar`i'
		bysort `id': gen `LIdepvar`i'' = `ldepvar`i''*`interaction'
	}
	* Get the model
	local model `"`model' (`ddepvar`i'' `ldepvar`i'' `LIdepvar`i'' `dvars' `lvars' `DIvars' `LIvars' `diffshockvariables' `lagshockvariables' `DIshockvariables' `LIshockvariables' `dummy' `dumIs')"'
	local i = `i' + 1
}
local maxdv = `i'-1							// need the max # of compositions
qui sureg `model'							// run the model and keep sample
qui keep if e(sample)
noi estsimp sureg `model'					// run Clarify
* ------------------------ Scalars and Setx ---------------------
qui setx mean								// set everything to means to start

qui su `lagshockvar', meanonly				// scalars for lagged shock
local sv = r(mean)
local vs = `sv' + `shock'
qui forv i = 2/5	{						// and the additional shocks
	if "`shockvar`i''" != "" {
		su `lagshockvar`i'', meanonly
		local sv`i' = r(mean)
		local vs`i' = `sv`i'' + `shock`i''
		setx `lagshockvar`i'' mean
		setx `diffshockvar`i'' 0
	}
}

qui setx `diffshockvar' 0					// set differenced shock to 0
qui setx `lagshockvar' mean					// set lag shock to mean

qui if "`interaction'" != ""	{			
	if "`intype'" == "on"	{
		setx `LIshockvar' `sv'				// set Interacted shock var & others
		setx `DIshockvar' 0
		forv i = 2/5	{
			if "`shockvar`i''" != ""	{
				setx `LIshockvar`i'' `sv`i''
				setx `DIshockvar`i'' 0
			}
		}
		foreach var of varlist `varlist'	{
			su `L`var'', meanonly
			local `var'_mean = r(mean)
			setx `LI`var'' (``var'_mean')
			setx `DI`var'' 0
		}
	}
	else	{
		setx `LIshockvar' 0
		setx `DIshockvar' 0
		forv i = 2/5	{
			if "`shockvar`i''" != ""	{
				setx `LIshockvar`i'' 0
				setx `Dishockvar`i'' 0
			}
		}
		foreach var of varlist `varlist'	{
			setx `LI`var'' 0
			setx `DI`var'' 0
		}
	}
}


qui forv i = 1/`maxdv' { 						// set lagged DVs to their means
	 su `mdepvar`i'', meanonly
	 scalar m`i' = r(mean)
	 setx `ldepvar`i'' m`i'
	 if "`interaction'" != ""	{				// and the lagged inter DV
	 	if "`intype'" == "on"	{
			setx `LIdepvar`i'' m`i'
		}
		else	{
			setx `LIdepvar`i'' 0
		}
	 }
}

qui setx (`dvars') 0						// set diff indep vars to 0
qui setx (`lvars') mean

local i 1
qui if "`dummy'" != "" {						// set our dummies, if they exist
	if "`dummyset'" != "" {
		foreach var in `dummy' {
			local m 1
			foreach k of numlist `dummyset' {
				if `m' == `i' {
					setx `var' `k'
				}
				local m = `m' + 1
			}
			local i = `i' + 1
		}
	}
	else {
		setx(`dummy') 0
	}
}

local i 1
qui if "`dummy'" != ""	{	// set interactive dummies to mean of shock, if they exist & if dummy = 1
	if "`interaction'" != ""	{
		foreach var in `dummy'	{
			if "`interaction'" != ""	{
				if "`intype'" == "on"	{
					local m 1
					foreach k of numlist `dummyset'	{
						if `m' == `i' & "`k'" == "1"	{
							setx(`dumI`var'') 1		// dummy interaction ON
						}
						local m = `m' + 1
					}
				}
				else	{
					setx(`dumI`var'') 0				// Dummy interaction OFF
				}
			}
			local i = `i' + 1
		}
	}
}


* ------------------------ Predict Values, t = 1 ----------------------------------
local preddv
qui forv i = 1/`maxdv' {
	tempvar td`i'log1
	local preddv `"`preddv' `td`i'log1' "'
}
qui simqi, ev genev(`preddv')				// grab our expected values

local denominator1
local i 1
qui foreach var in `dvs' {
	su `mdepvar`i''
	scalar z = r(mean)
	tempvar t`i'log1
	gen `t`i'log1' = z + `td`i'log1'
	su `t`i'log1', meanonly
	scalar m`i' = r(mean)				// these scalars become the new LDV
	local denominator1 `"`denominator1' + (exp(`t`i'log1'))"' // need this for below
	local i = `i' + 1
}
* ------------------------Loop Through Time-------------------------------------
di ""
nois _dots 0, title(Simulation in Progress) reps(`range')
qui forv i = 2/`range' {
	noi _dots `i' 0
	
	forv x = 1/`maxdv' {
		setx `ldepvar`x'' (m`x')					// set the new value of LDV
   	 	if "`interaction'" != ""	{				// and the lagged inter DV
   	 		if "`intype'" == "on"	{
   				setx `LIdepvar`x'' m`x'
   			}
   			else	{
   				setx `LIdepvar`x'' 0
   			}
   	 	}
	}	 
	
	* Set dummy variables:
	local op 1
	qui if "`dummy'" != "" {						// set our dummies, if they exist
		if "`dummyset'" != "" {
			foreach var in `dummy' {
				local m 1
				foreach k of numlist `dummyset' {
					if `m' == `op' {
						setx `var' `k'
					}
					local m = `m' + 1
				}
				local op = `op' + 1
			}
		}
		else {
			setx(`dummy') 0
		}
	}
	local op 1
	qui if "`dummy'" != ""	{	// set interactive dummies to mean of shock, if they exist & if dummy = 1
		if "`interaction'" != ""	{
			foreach var in `dummy'	{
				if "`interaction'" != ""	{
					if "`intype'" == "on"	{
						local m 1
						foreach k of numlist `dummyset'	{
							if `m' == `op' & "`k'" == "1"	{
								setx(`dumI`var'') 1		// dummy interaction ON
							}
							local m = `m' + 1
						}
					}
					else	{
						setx(`dumI`var'') 0				// Dummy interaction OFF
					}
				}
				local op = `op' + 1
			}
		}
	}
	
	qui setx `lagshockvar' mean	
	qui setx `diffshockvar' 0
	qui if "`interaction'" != ""	{			
		if "`intype'" == "on"	{
			setx `LIshockvar' `sv'			// set Interacted shock var & others
			setx `DIshockvar' 0
			forv z = 2/5	{
				if "`shockvar`z''" != ""	{
					setx `LIshockvar`z'' `sv`z''
					setx `DIshockvar`z'' 0
				}
			} 	
		}
		else	{
			setx `LIshockvar' 0			
			setx `DIshockvar' 0
			forv z = 2/5	{
				if "`shockvar`z''" != ""	{
					setx `LIshockvar`z'' 0
					setx `DIshockvar`z'' 0
				}
			}
		}
	}
	
		
	qui if `i' == `time' {					// we experience the shock at t
		setx `diffshockvar' (`shock')	// shock affects at time t only
		forv l = 2/5	{				// and the additional shocks if != ""
			if "`shockvar`l''" != "" {
				setx `diffshockvar`l'' (`shock`l'')
			}
		}
		if "`interaction'" != ""	{
			if "`intype'" == "on"	{
				setx `DIshockvar' (`shock')
				forv l = 2/5	{
					if "`shockvar`l''" != ""	{
						setx `LIshockvar`l'' (`shock`l'')
					}
				}	
			}
		}
	}
	else if `i' > `time' {
		setx `diffshockvar' 0			// diff shock back to 0
		setx `lagshockvar' (`vs')		// lag shock now at (mean + shock)
		forv  l = 2/5	{
			if "`shockvar`l''" != ""	{
				setx `diffshockvar`l'' 0
				setx `lagshockvar`l'' (`vs`l'')
			}
		}
		if "`interaction'" != ""	{
			if "`intype'" == "on"	{
				setx `DIshockvar' 0
				setx `LIshockvar' (`vs')
				forv l = 2/5	{
					if "`shockvar`l''" != ""	{
						setx `DIshockvar`l'' 0
						setx `LIshockvar`l'' (`vs`l'')
					}
				}
			}
		}
	}
	else {
		setx `diffshockvar' 0			// just to be sure
		forv  l = 2/5	{
			if "`shockvar`l''" != ""	{
				setx `diffshockvar`l'' 0
			}
		}
		if "`interaction'" != ""	{
			if "`intype'" == "on"	{
				setx `DIshockvar' 0
				forv l = 2/5	{
					if "`shockvar`l''" != ""	{
						setx `DIshockvar`l'' 0
					}
				}
			}
		}
	}
	
	qui setx (`dvars') 0				// just to be sure
	qui setx (`lvars') mean
	if "`interaction'" != ""	{
		foreach var of varlist `varlist'	{
			if "`intype'" == "on"	{
				setx `LI`var'' (``var'_mean')
				setx `DI`var'' 0
			}
			else	{
				setx `LI`var'' 0
				setx `DI`var'' 0
			}
		}
	}	
	
	local preddv
	qui forv x = 1/`maxdv' {
		tempvar td`x'log`i'
		local preddv `"`preddv' `td`x'log`i'' "'
	}
	simqi, ev genev(`preddv')			// get new predictions
	local denominator`i'
	qui forv x = 1/`maxdv' {
		tempvar t`x'log`i'
		gen `t`x'log`i'' = m`x' + `td`x'log`i''	// add them to old predictions
		su `t`x'log`i'', meanonly
		scalar m`x' = r(mean)
		local denominator`i' `"`denominator`i'' + (exp(`t`x'log`i''))"' // need this for below
	}
}

* keep only start-time, shocktime and endtime:
loc keeps
forv m = 1/`maxdv'	{
	local keeps `"`keeps' `t`m'log`basetime'' `t`m'log`time'' `t`m'log`endtime'' "'
}
keep `keeps'


* ------------------------Un-Transform----------------------------------------
* note that after untransforming compositions, starting means are subtracted and
* percentiles drawn.
local keepthese
local z = `maxdv' + 1
	qui forv m = 1/`maxdv' {
		tempvar var_pie`basetime'_`m' var_pie`time'_`m' var_pie`endtime'_`m'
		* Create compositions for start time:
		gen `var_pie`basetime'_`m'' = (exp(`t`m'log`basetime''))/(1  `denominator`basetime'')
		_pctile `var_pie`basetime'_`m'', p(50)		// grab midpoint
		gen var_pie_med_`basetime'_`m' = r(r1)
		
		* create composition for shocktime:
		gen `var_pie`time'_`m'' = ((exp(`t`m'log`time''))/(1  `denominator`time'')) - var_pie_med_`basetime'_`m'
		_pctile `var_pie`time'_`m'', p(`sigl',`sigu')
		gen var_pie_ll_sr_`m' = r(r1)
		gen var_pie_ul_sr_`m' = r(r2)
		
		* create composition for LR-time:
		gen `var_pie`endtime'_`m'' = (exp(`t`m'log`endtime''))/(1  `denominator`endtime'') - var_pie_med_`basetime'_`m'
		_pctile `var_pie`endtime'_`m'', p(`sigl',`sigu')
		gen var_pie_ll_lr_`m' = r(r1)
		gen var_pie_ul_lr_`m' = r(r2)
		local keepthese `"`keepthese' var_pie_ll_sr_`m' var_pie_ul_sr_`m' var_pie_ll_lr_`m' var_pie_ul_lr_`m' "'
	}				  
	tempvar var_pie`time'_`z' var_pie`endtime'_`z' var_pie`basetime'_`z'  // the un-transformation baseline
	gen `var_pie`basetime'_`z'' = 1/(1  `denominator`basetime'')
	_pctile `var_pie`basetime'_`z'', p(50)		// grab midpoint
	gen var_pie_med_`basetime'_`z' = r(r1)
	
	gen `var_pie`time'_`z'' = (1/(1  `denominator`time'')) - var_pie_med_`basetime'_`z'	
	_pctile `var_pie`time'_`z'', p(`sigl',`sigu')
	gen var_pie_ll_sr_`z' = r(r1)
	gen var_pie_ul_sr_`z' = r(r2)
	gen `var_pie`endtime'_`z'' = (1/(1  `denominator`endtime'')) - var_pie_med_`basetime'_`z'
	_pctile `var_pie`endtime'_`z'', p(`sigl',`sigu')
	gen var_pie_ll_lr_`z' = r(r1)
	gen var_pie_ul_lr_`z' = r(r2)
	
	local keepthese `"`keepthese' var_pie_ll_sr_`z' var_pie_ul_sr_`z' var_pie_ll_lr_`z' var_pie_ul_lr_`z' "'

keep `keepthese'
qui keep in 1

gen count = _n
reshape long var_pie_ll_sr_ var_pie_ul_sr_ var_pie_ll_lr_ var_pie_ul_lr_, j(sort) i(count)
gen sort_lr = sort - 0.2
gen sort_sr = sort + 0.2


* create midpoints:
gen mid_sr = (var_pie_ll_sr_ + var_pie_ul_sr_)/2
gen mid_lr = (var_pie_ll_lr_ + var_pie_ul_lr_)/2

/*
tempvar count 
qui gen `count' = _n

local reshapevar
local maxvars = `maxdv' + 1					// get max number of DVs
qui forv i = 1/`maxvars' {
	local reshapevar `"`reshapevar' var`i'_pie_ll_ var`i'_pie_ul_ "'
}

qui reshape long `reshapevar', j(time) i(`count')
qui drop `count'

qui forv i = 1/`maxvars' {
	gen mid`i' = (var`i'_pie_ll_ + var`i'_pie_ul_)/2
}
*/
di ""
if "`saving'" != "" {
	noi save `saving'.dta, replace
}
else	{
	noi save dynsimpie_results.dta, replace
}

end 										// end
* ----------------------------------------------------------------------------
